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Abstract. Multiconfiguration expansions frequently target valence correlation and 
correlation between valence electrons and the outermost core electrons. Correlation 
within the core is often neglected. A large orbital basis is needed to saturate both 
the valence and core-valence correlation effects. This in turn leads to huge numbers 
of CSFs, many of which are unimportant. To avoid the problems inherent to the 
use of a single common orthonormal orbital basis for all correlation effects in the 
MCHF method, we propose to optimize independent MGHF pair-correlation functions 
(PCFs), bringing their own orthonormal one-electron basis. Each PCF is generated 
by allowing single- and double- excitations from a multireference (MR) function. This 
computational scheme has the advantage of using targeted and optimally localized 
orbital sets for each PCF. These pair-correlation functions are coupled together and 
with each component of the MR space through a low dimension generalized eigenvalue 
problem. Nonorthogonal orbital sets being involved, the interaction and overlap 
matrices are built using biorthonormal transformation of the coupled basis sets followed 
by a counter-transformation of the PCF expansions. 

Applied to the ground state of beryllium, the new method gives total energies that 
are lower than the ones from traditional CAS-MCHF calculations using large orbital 
active sets. It is fair to say that we now have the possibility to account for, in a 
balanced way, correlation deep down in the atomic core in variational calculations. 
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1. Introduction 

There are several widely used methods for calculating many-electron atoms, but many- 
body effects still represent a real challenge. The many-body perturbation theory 
(MBPT) appears to be the most effective for atoms with one valence electron but its 
accuracy is not satisfactory for the atoms with more than one valence electron, mainly 
because of the poor convergence of the MBPT for the valence- valence correlations [1]. 
However, the core-valence correlations still can be effectively treated with MBPT. 
For this reason it was suggested to combine many-body perturbation theory for the 
core-valence correlations with the configuration interaction (CI) for valence-valence 
correlations within the (CI + MBPT) method [2], but this extension to atoms with 
more than three electrons in open shells met some difficulties. A method based on 
the so-called y^-^^ approximation has been proposed by Dzuba and Flambaum [3] 
including core-valence correlations by means of MBPT. The most recent extension of 
MBPT is the development of a configuration-interaction plus all-order method for atomic 
calculations [1] . This is a theoretical method combining the all-order approach currently 
used in precision calculations of properties of monovalent atoms with the configuration 
interaction approach that is applicable for many-electron systems. 

The coupled-cluster (CC) approach [5] is an interesting alternative, summing 
up all orders and taking into account pair correlations. Interesting computational 
developments have been proposed [H [7], not only in physics but also in quantum 
chemistry (for a complete review, see |8]). Successful applications do exist [91 fTOl [TT]. 
with some limitations for atoms with more complicated electron structures than alkali 
atoms. According to [2J, the most obvious shortcoming of the CC method is the neglect 
of three-particle correlations. Moreover, it treats the valence-valence and core-valence 
correlations at the same level of approximation while the former is much stronger than 
the latter. 

In quantum chemistry, variational complete active space self-consistent field 
(CASSCF) methods are quite successful for describing small and medium-size molecules 
but are not sufficient when external (dynamic) correlation must be included |12j . 
The latter are treated through second order perturbation theory using a single or 
multireference state [H]. Nonperturbative variational methods treat many-body 
effects in an accurate way for valence electrons, in both nonrelativistic and relativistic 
schemes. In this line, multiconfiguration Hartree-Fock (MCHF) and multiconfiguration 
Dirac-Hartree-Fock (MCDHF), often combined with CI methods p!5l [T6] have been 
widely used for accurate calculations of many-electron atomic properties, focusing on 
valence and core-valence correlation of large atomic systems. The accuracy of CI is 
limited by the incompleteness of the set of configurations used, if the one-electron 
orbital basis is complete. A practical limitation is the number of possible configuration 
state functions (CSFs) that becomes so large for a many-electron atom that one has 
to select only a small fraction of them [17j. This is often done by neglecting core 
excitations or only including a limited number of them [T9] . The CI method 
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is by definition (orbital) basis-dependent, while the multiconfiguration methods are 
not, the CI problem being iteratively coupled to the orbital optimization [201 IS]. As 
discussed in section HI the shape of the resulting orbitals strongly depends on the type of 
correlation introduced through the multiconfiguration expansion [22j. This is due to the 
properties of the variational principle applied for deriving the MCHF/MCDHF equations 
to be solved, with the consequence that a set of radial distributions resulting from a 
given correlation model/expansion could become inadequate, or at least incomplete, for 
another model or for an extension of the original one. Other problems encountered in 
variational multiconfiguration calculations are discussed in the present work, illustrating 
the difficulty of optimizing a single orthonormal orbital basis set from which CSF 
expansions would describe efficiently all correlation effects for a given physical state, and 
would produce reliable expectation values for any operator other than the Hamiltonian. 

To avoid these problems inherent to the use of a single common orthonormal orbital 
basis for all correlation effects in the MCHF method, we propose to first optimize 
independent MCHF pair correlation functions (PCFs), bringing their own orthonormal 
one-electron basis. These PCFs are then coupled to each other through a low dimension 
generalized eigenvalue problem. A pioneer and inspiring work was done by Froese Fischer 
and Saxena [23l [2l] who introduced the separated-pair MCHF approach. The originality 
of the present study lies in the way the one-electron nonorthogonalities are treated 
for setting up the interaction matrix between the reference and the PCF spaces. For 
each coupling matrix element, we adopt the biorthogonal orbital transformations and 
CI eigenvectors counter-transformations, as originally proposed by Malmqvist [25] and 
later extended to the spherical atomic symmetry by Olsen et al. |26] . 

The PCF interaction approach is tested, in the nonrelativistic approximation, on the 
ground state of beryllium, using a multireference description for the valence electrons. 
The results are compared with the complete active space (CAS) MCHF method. 



2. The variational method 



The MCHF method, that is extended and explored with respect to nonorthogonalites 
in the current work, can be derived from the variational principle. A state ^ is an 
eigenstate of the Hamiltonian H if and only if the energy functional 

is left unchanged for any infinitesimal variation in the state at the point \E'. For the 
ground state or states that are the lowest of it's symmetry, the variational method gives 
a minimum principle 

^° ^ (2) 
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E\^] being the upper bound to the exact ground state energy Eq. Using a superposition 
ansatz 

M 

\^>) = Y,cm (3) 

i=l 

for describing the model state in the subspace of the M basis states $j, the variational 
parameters {q} can be determined by solving the generalized eigenvalue problem 

He = ESc (4) 

with ifjj = {f^i\H\(^j) and S'jj = being, respectively, the Hamiltonian and 

overlap matrices. The total energies are found as roots of the secular equation 

det(H-^S) = 0. (5) 

The use of the eigenfunctions {\['('^)} as a model subspace satisfying 

(v]>W|v]>(0) = , {¥''>\H\^'^''>) = e,6ki (6) 

is useful for the description of excited states since they according to the Hylleraas- 
Undheim theorem [271 [28] satisfy the conditions 

Ek<e,, Vfc = l,...,M. (7) 

That is, approximate eigenvalues obtained by diagonalizing the Hamiltonian in a 
subspace can only be stabilized when the latter is enlarged. 

3. The nonrelativistie multieonfiguration Hartree-Foek method 

Starting from the nonrelativistie Hamiltonian for an A^-electron system 

TV p n N 

i=l L *-l i<j 

the MCHF approach determines an approximate wave function for the state labeled 
'jLS of the form 

M 

\^i^LS)) = Y,c,m^,LS)), (9) 

i=l 

where 7 represents the dominant configuration and any additional quantum numbers 
required for uniquely specifying the state being considered. The CSFs {$(7jL5')} are 
built from a basis of one-electron spin-orbitals 

4){nlmims) = -P{nl; r)YimXO,<^)Xm,, (10) 
r 

where the radial distributions {P{nl; r)} are to be determined. By applying the 
variational principle one obtains a set of integro-differential MCHF equations 

lf^ + '^[Z- Y{nl ; r)] - - e„,,„A P{nl ; r) = ^X{nl ; r) + ^ e„,„,P(n7 ; r)(ll) 
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for the unknown radial distributions [12]. The equations are coupled to each other 
through the direct Y and exchange X potentials and the Lagrange multipliers eniyi- 
The Lagrange multipliers force the radial orbitals to be orthonormal within the same / 
subspace. Under these conditions the configuration state functions are orthonormal 

m^,LS)m^,LS)) = 6,,. (12) 

The mixing coefficients {q} appearing in the expansion over CSFs also enter in the form 
of the potentials and can be determined by solving the CI problem 

He = Ec (13) 

for the current set of radial distributions. The MCHF and CI problems are solved iter- 
atively until self-consistency is reached for the radial distributions and for the selected 
Cl-eigenvector. 

The multiconfiguration method incorporates an extension allowing radial nonorthonor- 
mality, which sometimes can be used to advantage. Evaluation of matrix elements when 
orbitals are nonorthogonal is complex, in general. In order to keep the energy expres- 
sions manageable a number of restrictions are imposed: the orbitals within each CSF 
are mutually orthogonal, there are at most two subshells in $(7jLS') containing spec- 
tator electrons whose orbitals are nonorthogonal to orbitals in the interacting function 
^{'jjLS), if all the spectator electrons with nonorthogonal orbitals have the same / value 
then there are at most two such electrons in each of ^{jiLS) and ^{'jjES), the use of 
nonorthogonal orbitals is such that ($(7jLS')|$(7jL5')) = 6ij. Here the term spectator 
refers to those electrons not directly involved in the interaction [15]. Despite these re- 
strictions, nonorthogonal orbitals have successfully been applied in a number of cases, 
describing major correlation effects in a very compact way [22l 129] . 

4. Correlation and spatial location of orbitals 

The MCHF method results from the application of the variational principle, and the 
solution depends strongly on the energy functional or CSF expansion used to derive 
the MCHF equations. The direct influence of the CSF expansion on the shape of 
the resulting radial orbital basis can be exploited to target specific effects, but could 
also bring undesirable distortions of the wave function for the property of interest. To 
illustrate these effects let us look at the beryllium ls^2s^ ^S* ground state. The valence 
pair-correlation function (PCF) is described by considering all possible excitations of 
the valence electron pair 

\Ayy) = ai\ls^ 2s^ ^S) 

+ J2(^n\ls^2sns ^S)+ J2 ani,n'i'\ls^nln'l' ^S). (14) 

n nl,n'V 

The core-valence PCF is obtained by promoting one core electron together with one 
valence electron 

\Kcy) = P^\ls^2s^ ^S) + (3ni,n'i'\ls2snln'l' ^S). (15) 

nl,n'l' 
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Finally, the core-core PCF takes into account the correlation within the core 
\Kcc)=li\ls^2s^ 'S) 

'S)+ Yl lni,n'i'\2s^nln'l' 'S). 



(16) 




Each of the above sums represent correlation between a particular pair of electrons, 
namely (2s, 2s) for valence correlation, (Is, 2s) for core- valence, and (Is, Is) for core-core 
correlation. For this reason they are called pair-correlation functions. As far as termi- 
nology is concerned, pair-correlation should be strictly used for double-excitations [30] . 
However, single-excitations are also included in our PCFs. Somewhat arbitrary we 
choose to include |ls^2sns ^S) in the valence function and |ls2s^ns ^S) in the core- 
core function. Other choices are possible for the single-excitations. In the above case 
there is only one reference configuration, but applied to a multireference the PCFs cap- 
ture major correlation effects in a very effective way [151 l3T] . 

The dependence of the orbitals on the type of correlation is illustrated in Figure 1. Here 
one clearly sees the contraction of the correlation orbitals when going from a valence to 
a core-core correlation calculation. Unfortunately, correlation effects are not additive, 
and interference occurs through multiple excitations. As we have described above one 
can easily target specific correlation effects through the CSF expansion. However, the 
resulting set of radial distributions is localized in such a way that it becomes inadequate 
for representing other types of correlation. In our four-electron system it is for example 
not a meaningful to use core-core correlation orbitals for describing valence correlation 
and vice versa. The problem with the space localization of the correlation orbitals be- 
comes more pronounced for larger systems with many electron subshells. 

Multiconfiguration expansions frequently target valence correlation and correlation 
between valence electrons and the outermost core electrons. Correlation within the 
core is neglected. Based on the combined valence and core-valence expansion a single 
orbital basis is determined. While this approach gives orthonormal configuration state 
functions, and allows standard methods to be used for the construction of radial matrix 
elements, it has some drawbacks. A large orbital basis is needed to saturate both the 
valence and core-valence correlation effects. This in turn leads to huge numbers of CSFs, 
many of which are unimportant. Another problem is uneven convergence patterns as 
new layers of orbitals, dependent on the gain in variational energy, may be more or less 
contracted, affecting computed properties differently. Finally, there are difficulties to 
effectively incorporate correlation within the core. A solution to these problems, taking 
the nonorthogonal extension of the MCHF package all the way, would be to represent 
the wave function as an expansion over CSFs 



i=l 

where each CSF, depending on the correlation effect it describes, is built on an orbital set 
optimized for this effect. Since different orbital sets are involved, the CSFs are no longer 



M 




(17) 
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Figure 1. Contraction of tlie correlation orbitals when going from valence to core- 
valence and core-core correlation MCHF calculations of Be Is^ 2s^ ^S. The two thick 
lines correspond to the spectroscopic Is (no node) and 2s (one node) orbitals. Other 
lines represents the radial distributions of the correlation orbitals of the n = 4 active 
set. 



orthonormal, and the expansion coefficients are obtained by solving the generahzed 
eigenvalue problem. The Hamiltonian matrix elements coupling CSFs built on different 
and nonorthogonal orbital sets, could be evaluated using their Slater determinant 
expansions and the Lowdin's cofactor method [321 [33] . as adopted in many atomic 
structure packages |M1 |35l EH [37]. For CSF expansions, a nonorthogonal extension 
of the underlying Racah algebra [38] is available [39l HQ] but is not general enough. 
Moreover, both approaches quickly become prohibitive from a computational point of 
view. 



5. Pair-Correlation function interaction calculations 

A different way of handling nonorthogonalities is offered by the biorthonormal 
transformation technique originally introduced by Malmqvist [25]. The transformation 
allows general matrix elements between PCFs to be computed fast and efficiently. To 
exploit this we represent the wave function, not as an expansion over CSFs, but as an 
expansion over one or more reference CSFs |$[) and a number of known pair-correlation 
functions |Aj) built on separate and optimally localized orbital sets 

m = j2^^m) + j2~^^^. as) 

i=i j=i 

To avoid redundancies in the representation, reference CSFs are discarded from the 
PCFs. In the remaining sections we will use the notation |Aj) for PCFs including 
the reference CSFs and |Aj) for PCFs where the reference CSFs have been discarded 
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(weights of reference CSFs have been set to zero). The Hamiltonian and overlap matrix 
elements between reference CSFs and PCFs and between different PCFs are computed 
using the biorthonormal transformation. The expansion coefficients {q, c^} are obtained 
by solving the corresponding generalized eigenvalue problem. The problem of computing 
a total wave function including valence, core-valence and core-core correlation thus 
reduces to a series of separate pair-correlation problems together with a relatively small 
generalized eigenvalue problem. Practically, the calculations proceed in the following 
steps: 

(i) Perform an MCHF calculation for the reference expansion. 

(ii) Keep the orbitals from the first step fixed and perform separate MCHF calculations 
for the different PCFs including the reference. 

(iii) Remove the reference CSFs from the PCFs by setting the weights to zero. 

(iv) Loop over CSFs and PCFs, perform biorthogonal transformations and evaluate the 
Hamiltonian and overlap matrix elements. 

(v) Solve the generalized eigenvalue problem. 

This computational scheme has the advantage of using targeted and optimally localized 
orbital sets for each of the PCFs. It is possible to include core correlation in a tractable 
way. Also, it is possible to determine the contribution to computed properties from each 
of the PCFs. What is lost is some variational freedom in the coefficients when going 
from an expansion of CSFs to and expansion over reference CSFs and PCFs. Below we 
describe the biorthogonal transformation at heart of the method. 

6. Biorthogonal transformations 

In the above context, solving the generalized eigenvalue problem (jl]) of dimension {g+p) 
in the basis of the g reference CSFs and p PCFs 

|Ai),|A2),...|A,)} (19) 

requires the calculation of the Hamiltonian and Gram (overlap) matrix elements. In 
the present approach, freezing the orbitals optimized for the MCHF reference function 
guarantees the orthononormality of orbitals subsets involved in {f^\\H\Kj) and ($^|Aj) 
matrix elements, but one-electron nonorthogonalities definitely appear in off-diagonal 
matrix elements coupling different PCFs, i.e. (Aj|iJ|Aj). In the most general case 
invoking other optimization strategies, one-electron nonorthogonalities will show up in 
all but the diagonal matrix elements. The biorthonormal approach is then applied for 
evaluating such matrix elements. 

The biorthogonal transformation is explained in details in [26] focusing on the 
calculation of transition probabilities using nonorthogonal orbitals. This method has 
shown to be very efficient and useful, producing reliable atomic transition data [HI W2\ 
through the atomic structure package ATSP2K [13]. It has also been implemented 
in GRASP2K [H] for the calculation of transition amplitudes in the full relativistic 
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context [in]. The idea is simple: two orbital sets that are not orthonormal to each 
other are first transformed to become biorthonormal. For a coupling matrix element 
(Ai|iJ|Ar) built in their own orbital basis {0f } and {0f } that are not orthonormal 

(0f I0f ) = Slf , (20) 
linear transformation^ 

0^ = 0^C^^ ; 0^ = 0^0^^, (21) 

are found to transform the two original orbital sets into two new biorthonormal sets 

(0f I0f ) = • (22) 
The advantage of the biorthonormality property (!22l) of the transformed orbital sets 
is that the evaluation of any matrix elements can proceed as in the orthonormal case, 
as originally found by Moshinsky and Seligman |35]. There is an infinity of pairs of 
transformation matrices (C"^"^, C"^^) that produce biorthonormal basis sets. In our 
approach [26], the choice adopted is predicted by the restrictions on the configuration 
state function spaces used for and A,,. We require the transformation matrices to be 
upper triangular, a suitable choice for estimating the effect of the orbital transformation 
on the mixing coefficients {af} — > {af } / {af } {^f } giving the two representations 
of the |A/) and \Kr) functions in both original and transformed (biorthonormal) basis 
sets 

i i 



Malmqvist |25] has indeed shown that an upper-triangular orbital transformation 
matrix can be expressed as a finite sequence of single-orbital transformations, each 
expressing the new orbitals as a sum involving no higher-numbered orbital. Each such 
transformation step on the Cl-expansion array is the same as the effect of a one-electron 
operator with de-excitations only. To avoid symmetry-breaking intermediates in the 
recursive transformation, some atomic symmetry refinement was made in [26] to express 



the transformation operator acting in the CSF space in the following suitable form 

^2(2i+l) \ 

Ei^'^b""' (25) 



N=0 



with 



'2(2«+l) 



(26) 



Nk is the occupation number of the (n/)fc-subshell and t is the matrix defining the single 
orbital replacements sequence that can be calculated from the UL decomposition of 

X in our matrix notation, the orbitals {(j)i} are collected in row vectors 4>. 
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(S^^) ^ . The knowledge of the action of the excitation operators in the CSF spaces for 
both |A/) and \Kr) 

5^ (^nnm.mAnlm,^. |$.) = j) (27) 

\ minis J j 

is then enough to perform the countertransformations of the corresponding CI vectors. 
Starting from the coupled tensorial second quantized form of the single-particle 
Hamiltonian operator 

H = -Y, \/2(2/ + 1) ^ (2^^ 

I n' ,n 

one realizes [26] that Aij = — a/2(2Z + 1) times the coefficient of the In'i,ni integral found 
in Hij (see also [IS])- 

There is an important built-in constraint in the algorithm: the wave function 
expansion spaces for both left and right should be "closed under de-excitation" to 
allow this class of transformation. Restricted active space (RAS) and complete active 
space (CAS) expansions |37] satisfy this property, i.e. removing one electron from a 
subshell nl and placing it in any subshell n'l of the same spatial symmetry [n' < n) 
generates a CSF that appears in the original configuration expansion. Other types of 
expansion can easily be complemented with CSFs to satisfy the closure condition. 



7. Code implementation 

The nonorthogonal pair-correlation method is implemented in a code module extending 
the ATSP2K package [13]. The module contains routines for checking closure under 
de-excitation (Iscud), performing biorthogonal transformations (biotrans), evaluating 
Hamiltonian and overlap matrix elements (biomatrix), and solving the generalized 
eigenvalue problem [18]. Having determined the reference and the pair-correlation 
functions, removing the reference CSFs from the latter, the calculation proceeds as 
follows: 

loop over reference and pair-correlation functions 

call Iscud 

call biotrans 

call biomatrix 
end loop 

call generalized eigenvalue solver 

To perform the biorthogonal transformation only one-electron coupling coefficients are 
needed, and thus this part is computationally fast. The subsequent evaluation of the 
Hamiltonian and overlap matrix elements utilizes standard Racah algebra techniques in 
the biorthonormal basis. Since there are two orbital sets being biorthogonal, rather than 
just one set, normal symmetry relations for radial integrals do not hold and integration 
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routines need to be re-designed. The correctness of the implemented method was checked 
by comparing Hamiltonian and overlap matrix elements with the corresponding ones 
calculated using Slater determinant algebra in nonorthogonal bases [321 |33]. Due to 
the complexity of the latter method this could only be done for two-, three- and four- 
electron systems. The correctness of the codes can also be inferred from calculations 
using artificially rotated orbital sets, where invariances in the representations are used. 
Tests for different types of systems showed that the biorthonormal transformation is 
numerically robust, with negligible loss of accuracy for the matrix elements. 

8. Exploring computational strategies 

The beryllium atom may be thought of as the "benchmark" atom. It is the first neutral 
system in which the three correlation effects, i.e. valence, core-valence and core corre- 
lation, appear and during the years much work has been devoted to understand and 
accurately describe these effects. Although the goal of the present work is not to get the 
lowest absolute total energy of the beryllium ground state, it is appropriate to highlight 
some important contributions with regard to chronology. Byron and Joachain |49j decou- 
pled the original four-electron problem into a series of helium-like equations describing 
pair correlation between electrons. A combined configuration-interaction- Hylleraas-type 
wave function study was performed by Sims and Hagstrom [50]. A rather complete cor- 
relation study of Be was presented by Froese Fischer and Saxena [23] who introduced 
the separated-pair MCHF approach (this latter work is undoubtedly the inspiration 
source of the present approach). Elaborate and highly accurate Be ground state to- 
tal energies were predicted by Bunge [5T|. An iterative numeric procedure [52] was 
used to obtain pair functions applied to two-electron systems. Numerical many-body 
perturbation calculations on Be-like systems were proposed by Salomonson et al. [53] 
adopting a multiconfigurational model space. A unified approach combining the multi- 
configuration Hartree-Fock method and many-body perturbation theory was attempted 
by Morrison [SU [311 ISS] . The coupled-cluster single- and double-excitation equations 
were solved numerically by Salomonson and Oster [56]. After the beryllium atom was 
"revisited" by a number of groups [57] , the nonrelativistic total energy of the Be ground 
state was estimated by Lindroth et al. [58j. Accurate values for the nonrelativistic 
total energy of Be ground state have been obtained using the full-core plus correlation 
method [59]. A new correlation study by Froese Fischer [17] was presented in the MCHF 
scheme, extending the study of n-expansion methods to four-electron systems. Energy 
levels and oscillator strengths were calculated by Weiss [60] adopting a multireference 
superposition-of-configurations approach for describing core-correlation effects. Finite- 
element MCHF and GTO basis set expansions coupled with MRSDCI were used to 
estimate the beryllium electron affinity by Olsen et al. [61] . Although difficulties appear 
in computing matrix elements in the Hylleraas basis [62l [63] , Hylleraas configuration- 
interaction calculations were performed on the nonrelativistic ground-state energy of 
the Be atom [6l]. The most accurate energies for the ground state of the beryllium- 
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like atoms have been obtained using the exponentially correlated Gaussian basis sets 
[65| |66] . The most recent work [66] reports the up-to-date lowest infinite mass nonrela- 
tivistic total energy of —14.667 356 486 E^^. 

Below we explore different computational strategies for ls^2s^ in beryllium 
based on separately optimized pair-correlation functions. The results are compared 
with calculations utilizing a single ortho normal orbital set. Whereas many accurate 
computational schemes can only be applied to small systems, the current method is 
directly generalizable to more complex systems and cases for which it currently is not 
possible to saturate the orbital basis (see section [9]). 

Monoreference ls^2s^ - {A x 4) approach 

We start by generating the three separate PCFs: the valence, the core-valence, and the 
core-core 

\Avv) = ai\ls^ 2s^ ^S) 

+ ^a„|ls^ 2s ns ^S) + a^i^n'i'lls'^ nl n'l' ^S) (29) 

n nl,n'l' 

\Acv)= Pills' 2s' 'S) 

+ Yl Pni,n'i'\ls2snl n'l' ^S) (30) 

nl,n'l' 

\Acc)= nils' 2s' 'S) 

+ ^7„|ls 2s' ns ^S) + J2 lni,n'i'\2s' nl n'l' 'S) (31) 

n nl,n'V 

The Is and 2s orbitals are taken from an initial HF calculation and kept frozen. All 
correlation orbitals are variational resulting in three sets of orbitals. The energies and 
number of CSFs as function of the largest principal quantum number are displayed in 
Tabled] for each of the PCFs. The table shows, as expected, larger correlation energies 
for valence and core excitations, relatively to core- valence. 

Denoting the PCFs in which the reference CSF |ls^2s^ ^S) has been removed by 
setting the corresponding expansion coefficient to zero by, respectively, |Ayy), \kcv)) 
and I Ace), we write the wave function as 

1^) = ci,22,2|ls2 2s' ^S) + cvv\~^vv) + ccv\~kcv) + ccc\~kcc)- (32) 

The expansion coefficients and the total energies are obtained by solving the 
corresponding 4x4 generalized eigenvalue problem, see Table El The total energies 
are compared with the SD-MCHF results obtained from the "traditional" approach, i.e. 
- i) generating the CSF expansion by single- and double-excitations from the reference 
state, and - ii) varying a common set of the orbitals representing the three correlation 
effects. The convergence with respect to the increasing orbital set is initially much 
faster using the nonorthogonal pair-correlation approach, but as the orbital basis start 
to saturate the different correlation effects, the traditional method gives a lower total 
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Table 1. Energies together with the number of CSFs for the valence (VV), core- valence 
(CV), and corc-corc (CC) PCFs. The (Is, 2s) orbitals arc from HF. Correlation orbitals 
for each of the PCFs are determined in separate MCHF calculations. 



n < 




NcSF 




NcSF 


Ecc 


^CSF 


HF 


-14.573 023 17 


1 


-14.573 023 17 


1 


-14.573 023 17 


1 


2 


-14.616 062 66 


2 


-14.574 237 22 


2 


-14.595 158 32 


2 


3 


-14.618 619 53 


7 


-14.578 062 34 


6 


-14.612 266 39 


7 


4 


-14.618 990 08 


16 


-14.578 781 57 


19 


-14.614 398 07 


16 


5 


-14.619 083 21 


30 


-14.578 960 23 


40 


-14.615 041 50 


30 


6 


-14.619 121 18 


50 


-14.579 021 81 


72 


-14.615 315 90 


50 


7 


-14.619 138 74 


77 


-14.579 050 02 


117 


-14.615 455 15 


77 


8 


-14.619 148 03 


112 


-14.579 064 66 


177 


-14.615 530 36 


112 


9 


-14.619 153 81 


156 


-14.579 071 26 


254 


-14.615 580 19 


156 


10 


-14.()19 157 24 


210 


-14..")79 ()7() 87 


:-!."•)() 


-14.015 (ill 88 


210 



Table 2. Solution of the (4 x 4) generalized eigenvalue problem. The energies are 



couii)ar(xl witli tl 


w SD-AICIIF rcsullH 


Ijascd on a 




ortlioiioruial oil)! 


al sol. 


n < 




cvv 


ccv 




ccc 




-£^4x4 


EsD-MCHF 


2 


0.955006 


0.294901 


1.1855[- 


2] 


2.9271[ 


-2] 


-14.636 852 03 


-14.616 852 26 


3 


0.959405 


0.278621 


2.1230[- 


2] 


3.8235[ 


-2] 


-14.658 281 46 


-14.651 174 28 


4 


0.959708 


0.277380 


2.2430[- 


2] 


3.8948[ 


-2] 


-14.661 298 22 


-14.657 844 72 


5 


0.959832 


0.276919 


2.2651[- 


2] 


3.9046[ 


-2] 


-14.662 171 17 


-14.660 930 03 


6 


0.959920 


0.276607 


2.2714[- 


2] 


3.9061[ 


-2] 


-14.662 527 85 


-14.661 861 16 


7 


0.959978 


0.276406 


2.2737[- 


2] 


3.9062[ 


-2] 


-14.662 704 55 


-14.662 413 74 


8 


0.960006 


0.276305 


2.2746[- 


2] 


3.9059[ 


-2] 


-14.662 799 45 


-14.662 742 11 


9 


0.960024 


0.276246 


2.2750[- 


2] 


3.9057[ 


-2] 


-14.662 858 90 


-14.662 908 34 


10 


0.960035 


0.276206 


2.2752[- 


2] 


3.9055[ 


-2] 


-14.662 897 85 


-14.663 013 15 



energy. Looking at the results one should bear in mind that, for a given active set 
specified by the largest principal quantum number n, the number of variational orbitals 
is larger by a factor three in the pair-correlation approach, while the size of the PCF 
expansions are smaller than the traditional SD-MCHF expansion by roughly the same 
factor. 

Multireference Is^ {2s, 2p}^ - {5 x 5) approach 

In this section we investigate a multireference calculation based on the Layzer's complex 
Is^ {2s^ + 2p^} ^S. The CSF list is generated by making simple and double excitations 
on this complex to an active set of orbitals. The CSFs are then arranged into the three 
lists, one for each correlation type, 

|A^^) = ai\ls^ 2s^ ^S) + aalls^ 2p2 ^S) 

+ an I Is^ 2s ns ^S) -\- am \ Is^ 2p rap ^S) 
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+ ^ ani,n'i'\ls^ nln'l' ^S) (33) 

nl,n'V 

\kcv) = Alls' 2s^ ^S) + p2\ls^ 2p^ ^S) 
+ (^ni,n'i'\ls 2snl n'l' 'S) 

nl,n'l' 

+ E /5;.,nT|ls2pn/n'r (34) 

nl,n'l' 

\Acc)= nils' 2s' 'S)+^2\ls' 2p''S) 

+ J2ln\ls 2s' US ^S) + ^7^113 2/ mp ^S) 

n m 

+ IniM'^s' nl n'l' 'S) + Y <WI'\2P' n'l' 'S). (35) 

nl,n'l' nl,n'l' 

This kind of formal writing shows us that the CSF hst corresponding to the valence 
correlation is rigorously the same as the one in the single reference case. However 
the content of the two other lists change, they now include simple, double, triple and 
quadruple excitations in comparison with the fundamental configuration ls^2s^ ^S. 
Table [3] contains the results of the three independent pair-correlation MCHF calculations 
in which the Is, 2s, 2p orbitals are frozen from the calculation of the multireference. 
The remaining correlation orbitals are completely variational in all the calculations. 

Table 3. Energies together with the number of CSFs for the multireference 
ls2 {2s,2p}2 is* valence (VV), core-valence (CV), and core-core (CC) PCFs. The 
Is , 2s, 2p orbitals are from the multireference MCHF calculation. Correlation orbitals 
(except 2p) for each of the PCFs are determined in separate MCHF calculations. 



n < 




NcSF 




NcSF 


Ecc 


NcSF 


MR 


-14.616 845 32 


2 


-14.616 845 32 


2 


-14.616 845 32 


2 


3 


-14.618 905 91 


7 


-14.619 426 25 


12 


-14.654 311 44 


22 


4 


-14.619 124 07 


16 


-14.621 265 76 


40 


-14.658 148 75 


58 


5 


-14.619 200 57 


30 


-14.621 676 31 


93 


-14.659 044 75 


119 


6 


-14.619 234 08 


50 


-14.621 766 58 


177 


-14.659 373 26 


211 


7 


-14.619 251 00 


77 


-14.621 801 87 


298 


-14.659 517 62 


340 


8 


-14.619 260 03 


112 


-14.621 819 73 


462 


-14.659 609 71 


512 


9 


-14.619 265 73 


156 


-14.621 828 19 


675 


-14.659 662 96 


733 


10 


-14.619 268 69 


210 


-14.621 833 99 


943 


-14.659 696 30 


1009 



Again, using the notation lAyy), lAcv-), and \ Acc) for the PCFs where the weights 
of the reference CSFs have been set to zero, the wave function is written 

1^) = ci,22,2|ls2 2s2 ^S) + Ci,22p2|ls2 2p2 ^S) 

+ cvv\Avv) + ccv\Acv) + ccc\Acc)- (36) 

The expansion coefficients and the total energies are obtained by solving the 5x5 
generalized eigenvalue problem. The coefficients and the energies are reported in Table H] 
as functions of the largest principal quantum number in the expansions. The energies are 
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compared with values from traditional SD-MR-MCHF and CAS-MCHF calculations. 
In the 4x4 approach |ls^2p^ ^S) entered the lAyy) valence PCF, which obtained 
a comparatively large weight. Separating out |ls^2p^ ^S) increases the variational 
freedom through the expansion coefficient 01^22^2 and the valence PCF becomes just a 
small correction. Looking at the energies we see a dramatic improvement. The present 
5x5 energy for n = 10 is comparable to the n = 10 CAS-MCHF energy based on an 
expansion of more than 650 000 CSFs. The effort for the CAS-MCHF calculation is 
huge, and the case has to run for days on a cluster. In contrast, generating the pair- 
correlation functions, constructing the Hamiltonian and overlap matrices, and solving 
the 5x5 generalized eigenvalue problem is very fast and is easily done on a PC. 



Table 4. Solution of the (5 x 5) generalized eigenvalue problem. The energies are 
compared with SD-MR-MCHF and CAS-MCHF results based on a single orthonormal 
orbital set. 



n < 


Cls22s2 




Cvv 


Ccv 




Ccc 




3 


0.952146 


0.299633 


4.2579[-2] 


1.7789[- 


-2] 


3.8845[ 


-2] 


4 


0.952674 


0.297365 


4.3498[-2] 


2.1246[- 


'2] 


4.0538[ 


-2] 


5 


0.952773 


0.296950 


4.3639[-2] 


2.1823[- 


-2] 


4.0768[ 


-2] 


6 


0.952832 


0.296749 


4.3642[-2] 


2.1922[- 


'2] 


4.0805[ 


-2] 


7 


0.952864 


0.296644 


4.3629[-2] 


2.1956[- 


'2] 


4.0817[ 


-2] 


8 


0.952885 


0.296577 


4.3615[-2] 


2.1968[- 


'2] 


4.0816[ 


-2] 


9 


0.952897 


0.296543 


4.3603[-2] 


2.1972[- 


'2] 


4.0814[ 


-2] 


10 


0.952902 


0.296527 


4.3597[-2] 


2.1973[- 


-2] 


4.0813[ 


-2] 


n < 


-^5x5 


Emr-sd-mchf 


EcAS-MCHF 










3 


-14.658 887 01 


-14.654 399 79 


-14.654 414 59 










4 


-14.664 774 35 


-14.661 865 68 


-14.661 403 17 










5 


-14.666 173 94 


-14.664 721 18 


-14.664 839 93 










6 


-14.666 628 18 


-14.665 938 59 


-14.666 067 32 










7 

8 


-14.666 826 48 
-14.666 945 78 


-14.666 407 90 
-14.666 722 13 


-14.666 541 14 
-14.666 857 41 










9 


-14.667 013 54 


-14.666 876 64 


-14.667 012 75 










10 


-14.667 056 14 


-14.666 975 55 


-14.667 114 52 











Multireference Is^ {2s, 2p, 3s, 3p, 3(i}^ ^S* - (8 x 8) approach 

Here we extended the multireference to the n = 3 complex Is^ {2s^ -|- 2p^ + 3s^ -|- 3p^ + 
?)d?} ^S. As in the previous section we treat each component of the multireference 
as a subspace of the interaction matrix, increasing the dimension of the generalized 
eigenvalue problem to 8. The results of the separate pair-correlation MCHF calculations 
are displayed in Table [51 Even if the size of the largest expansion has increased, it is 
still a very small problem. Comparing with Table [3], we see that it is the core- valence 
and core correlation PCFs that now give lower energies. 
In this case the wave function has the form 

|^)= cmAIs^uI^^S) 

nl,n<3 
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Table 5. Energies together with the number of CSFs for the multireference 
Is^ {2s, 2p, 3s, 3p, 3d}^ valence (VV), core-valence (CV), and core-core (CC) 
PCFs. The n/, [n < 3) orbitals are from the multireference MCHF calculation. 
Remaining correlation orbitals for each of the PCFs are determined in separate MCHF 
calculations. 



n < 




NcSF 




NcSF 


Ecc 


NcSF 


MR 


-14.618 914 67 


5 


-14.618 914 67 


5 


-14.618 914 67 


5 


4 


-14.619 126 28 


16 


-14.622 479 34 


70 


-14.656 953 43 


153 


5 


-14.619 201 47 


30 


-14.623 667 44 


192 


-14.660 544 82 


352 


6 


-14.619 234 70 


50 


-14.623 918 77 


403 


-14.661 285 99 


667 


7 


-14.619 251 52 


77 


-14.623 977 10 


721 


-14.661 532 51 


1122 


8 


-14.619 260 43 


112 


-14.623 999 33 


1164 


-14.661 642 87 


1743 


9 


-14.619 266 16 


156 


-14.624 009 07 


1750 


-14.661 702 00 


2555 


10 


-14.619 269 59 


210 


-14.624 015 40 


2497 


-14.661 737 42 


3583 



+ cvv\-^vv) + ccv\-^cv) + ccc\-^cc)- (37) 

The expansion coefficients and the total energies are obtained by solving the 8x8 
generalized eigenvalue problem. The coefficients and the energies are reported in Table [6] 
as functions of the largest principal quantum number n in the expansions. The energies 
are compared with values from traditional CAS-MCHF calculations, which give the 
lowest possible energies that can be obtained from a single orthogonal orbital set. Even 
for n = 10, that is a basis consisting of 53 orbitals in the orthogonal case, the pair- 
correlation approach gives a lower total energy. The convergence with respect to the 
largest principal quantum number n is graphically displayed in figure O Here one clearly 
sees the slow saturation rate of the energy for the CAS-MCHF calculation. For larger 
systems with more subshells, the difference in saturation rate between calculations using 
nonorthogonal PCFs and traditional calculations built on a single orbital set will be 
much greater. Note that although the PCF interaction total energy is lower than the 
CAS-MCHF one, it is still 2.10~'^E-^ above the Stanke et a/.'s result [66]. Remembering 
that both the PCF- and CAS-MCHF expansions are /-truncated (/ < 10), this is 
expected due to the slow convergence rate {Ei — Ei^i 0{l + 1/2)^^) with respect 
to I [67|. 
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Table 6. Solution of the (8 x 8) generalized eigenvalue problem. The energies are 
compared with CAS-MCHF results based on a single orthonormal orbital set. 



n \ 




Cls2 2p2 






Cls23p2 


Cls23d2 


4 


u.yooiUo 


u.zyoyoD 


— 4.UDDo[ 




O.OZD4[— oj 


1 71 /1 r 01 
— 1. / i4y [— zj 





u.yooziu 


u.zyoz ( 


A ncncr 
— 4.UoUo[ 


01 


£).iDzy[— oJ 


— i.DyyD[— zj 







u.zyouUD 


— 4.U0U / [ 




o.ii /o[ — Oj 


— i.DoDy[ — ZJ 


7 
/ 


u.yooozz 


u.zy4oDy 


A n7CQr 
— 4.U( oo[ 


01 


iJ.UoU ( [—o\ 


— i.Doiy [— ZJ 


o 
o 


U.y0oo4z 


U.zy4oU4 


— 4.U /oy [ 


01 


0.UDZD[— oJ 


— i.D/oy[— ZJ 


9 


0.953353 


0.294771 


-4.0763[ 


-2] 


5.0516[-3] 


-1.6772[-2] 


10 


0.953362 


0.294743 


-4.0758[ 


-2] 


5.0409[-3] 


-1.6759[-2] 


n < 


cvv 


ccv 


ccc 




Esxs 


EcAS-MCHF 


4 


8.5613 [-3] 


2.0715 [-2] 


3.8931[- 


2] 


-14.660 679 48 


-14.661 403 17 


5 


9.6715[-3] 


2.1911[-2] 


4.0595[- 


2] 


-14.665 553 46 


-14.664 839 93 


6 


9.8881 [-3] 


2.2103[-2] 


4.0808[- 


2] 


-14.666 582 83 


-14.666 067 32 


7 


9.9702[-3] 


2.2135[-2] 


4.0844[- 


2] 


-14.666 905 87 


-14.666 541 14 


8 


9.9847[-3] 


2.2148[-2] 


4.0849[- 


2] 


-14.667 047 86 


-14.666 857 41 


9 


9.9822[-3] 


2.2151[-2] 


4.0848[- 


2] 


-14.667 122 76 


-14.667 012 75 


10 


9.9793[-3] 


2.2151[-2] 


4.0846[- 


2] 


-14.667 168 08 


-14.667 114 52 



-14.656 



-14.658 



-14.66 



-14.662 



-14.664 



-14.666 



-14.668 



\ MR Is^ {2s,2p,3s,3p,3d) S - 8x8 approach 1- 

\ CAS-MCHF — X- 




_1_ 



_±_ 



5 6 7 8 
maximum n value 



10 



11 



Figure 2. Energies of the Is^ {2s, 2p, 3s, 3p, 3^}^ multireference plus VV, CV, CC 
pair-correlation and the CAS-MCHF expansions as functions the principal quantum 
number n. 
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9. Conclusion and perspectives 

Variational calculations based on "explicitly correlated Gaussians" are quite impressive 
but unfortunately scale as A^!, where N is the number of electrons. As mentioned by 
Stanke et al. [66], calculations with fully correlated basis functions would require huge 
amounts of computer time for systems with more than six electrons. For the latter, 
the numerical MCHF method is widely used, generally limited to the optimization of 
a single orthonormal set of orbitals. However, the relaxation of one-electron orbitals 
and orthogonality constraints has been shown to be beneficial in many applications 
[68t [69t [29| 170] . The present work demonstrates the effectiveness of relaxing one-electron 
orthogonality constraints in the MCHF orbital optimization of different PCFs, yielding 
a new method for treating correlation. Applied to the ground state of beryllium the new 
method gives total energies that are lower than the ones from traditional CAS-MCHF 
calculations using large single orbital sets with principal quantum numbers up to 10. 
Whereas many accurate computational schemes can only be applied to small systems, 
the current method is directly generalizable to more complex systems and cases for which 
it currently is not possible to saturate a single orbital basis for describing different types 
of correlation contributing to the total energy, or different type of operators. It is fair 
to say that we now have the possibility to account for, in a balanced way, correlation 
deep down in the atomic core in variational calculations. 

In this study we have only looked at a few aspects of separately optimized PCFs. 
Following Froese Fischer and Saxena [231 El] one could further refine the method and 
define PCFs for each two-electron coupling. One would then, for example, differentiate 
between 

\Kcv) = ai\ls^2s^ ^S) + ^ ani,n'i'\ls2s ^S{nln'l' ^S) ^S) (38) 

nl,n'l' 

and 

\Acv) = Pi\ls'2s^ 'S) + J2 Pniyi'\ls2s ^S{nln'l' ^S) 'S). (39) 

nl,n'l' 

It would also be possible to introduce some single-electron excitation functions for 
describing core-opening effects crucial for hyperfine structure and other one-electron 
operator quantities. The method with separately optimized pair-correlation functions 
can be extended to a spectrum, describing several states at the same time. The atomic 
state function is then given by the superposition 

M / gr Pr \ 

i^) = E E<i^^) + Es^i^^) ■ (40) 

r=l \i=l j=l / 

Energies and expansion coefficients of the states are obtained by selected solutions of 
the eigenvalue problem. 

A lot of effort has been put [18] in the general adaptation of ATSP2K [43] for 
allowing separately optimized PCFs, and it is now possible to compute the expectation 
of any one- and two-electron operator, including the Breit-Pauli corrections, in the 
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new wave function representation. The biorthogonal transformation for handling 
nonorthogonahties is apphcable also in the fully relativistic case. This method with 
separately optimized PCFs is currently implemented [IE] in GRAPS2K [44j. 
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